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We extend the solid-state nudged elastic band method to handle a non-conserved order param¬ 
eter - in particular, magnetization, that couples to volume and leads to many observed effects in 
magnetic systems. We apply this formalism to the well-studied magneto-volume collapse during the 
pressure-induced transformation in iron - from ferromagnetic body-centered cubic (bee) austenite 
to hexagonal close-packed (hep) martensite. We find a bcc-hcp equilibrium coexistence pressure 
of 8.4 GPa, with the transition-state enthalpy of 156 meV/Fe at this pressure. A discontinuity in 
magnetization and coherent stress occurs at the transition state, which has a form of a cusp on the 
potential-energy surface (yet all the atomic and cell degrees of freedom are continuous); the calcu¬ 
lated pressure jump of 25 GPa is related to the observed 25 GPa spread in measured coexistence 
pressures arising from martensitic and coherency stresses in samples. Our results agree with exper¬ 
iments, but necessarily differ from those arising from drag and restricted parametrization methods 
having improperly constrained or uncontrolled degrees of freedom. 


I. INTRODUCTION 

Magneto-structural transformations are quite common 
in magnetic systems [T]. Their generic feature is a rapid 
change of magnetization and density, referred to as a 
magneto-volume collapse [2]. While the nudged elastic 
band (NEB) [3] and the solid-state nudged elastic band 
(SSNEB) |4 methods correctly account for all atomic 
degrees of freedom (DoF), they expect the transition 
state (TS) region in the form of a saddle. The oft-used 
drag methods can miss the correct TS (due to a possi¬ 
ble discontinuity in some of the DoF), and a restricted 
parametrization can overlook it (due to a constrained 
search space). Here we generalize the TS search algo¬ 
rithm within the SSNEB method to make it suitable 
for the TS of various shapes (not only saddles, but also 
cusps), and apply it to iron under pressure, which dis¬ 
plays a magneto-volume collapse of 11% at the TS due 
to the loss of magnetization pQ, a non-conserved order 
parameter that does not have to be continuous, unlike 
the atomic or cellular DoF. 

Iron is the most stable element produced by nuclear 
reactions, the most abundant element in the Earth core, 
and is also very common in extraterrestrial meteorites 
[5]. Metallic iron is known for its magnetism and tech¬ 
nological impact. Deeper understanding of its properties 
under pressure affects metallurgy, materials science and 
engineering, geophysics, and planetary sciences [6]. 

Due to its importance, the iron bcc-hcp ( as ) tran¬ 
sition and its mechanism have long been studied [6H2Z]. 
In spite of iron’s structural simplicity under pressure, the 
minimum-energy pathway (MEP) and the TS remained 
unresolved, in particular, because of the improper han¬ 
dling of magnetic DoF and expectation of saddle point 
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behavior near the TS. We use the bcc-hcp transition in 
iron as a prototype for the quantitative description of 
magnetic systems via the SSNEB method [4] and its mod¬ 
ifications [28]. For a magneto-structural transformation, 
we show that the TS is a cusp on the potential energy 
surface (PES) and is not a single electronic state, but a 
duality of a non-magnetic and magnetic states with the 
same atomic and cell coordinates, the same enthalpy, but 
different electronic structure. As a result of the coherency 
stress between those two electronic states, a pressure dis¬ 
continuity at the TS exists, which can be related to the 
experimental scatter in the coexistence pressures of mag¬ 
netic and non-magnetic phases. 

As discussed below, the simplest bcc-hcp transition in 
iron can be described in a 2-atom cell by one atomic 
and three cell (4 total) degrees of freedom, with the rest 
being constant during the whole transformation due to 
symmetry (Fig. [l]). The atomic DoF in a 2-atom cell can 
be completely described by a continuous shuffle s (which 
linearly changes from 0 for bee to 1 for hep); with direct 
coordinates of atom 1 fixed at (0; 0; 0) and atom 2 at 

(It 1 + fl; It 1 + fl; §)• 



FIG. 1: (color online) Left: Smallest (2-atom) bee unit cell 
that permits a bcc—)>hcp transform using a shuffle-shear con¬ 
cept. Right: bcc—)>hcp via shuffle-shear; bcc ([110] projection) 
is shown by black atoms and solid lines, and hep ([0001] pro¬ 
jection) is shown as blue open circles and dashed lines, with 
ABAB stacking of corner A and central B atoms. Vector s 
(red) gives the shuffle direction from bcc to hep. The shear 
angle d changes from 70° in bcc to 60° in hep. 
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FIG. 2: (color online) For 2-atom cell, enthalpy H [meV per 
atom] relative to bcc Fe, magnetization M [jib per atom], 
and volume V [A 3 per unit cell] vs. shuffle from SSNEB 
at 3 hydrostatic pressures: P=0 (black), 10 (red), and 20 
GPa (blue). Higher-enthalpy states with intermediate M are 
shown as filled shapes. Dotted vertical lines and dashes are 
at discontinuities in M. 
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FIG. 3: (color online) Lattice parameters a, b = 2a sin($/2), 
and c (A) versus shuffle for 3 values of external pressure 
(GPa): P —0 (black), 10 (red), and 20 (blue). Trajectory 
of the TS is given by dashed line. 

The 3 independent cell DoF for an orthorhombic unit 
cell are the lengths of the 3 mutually orthogonal vectors 
ai+a 2 , a 2 —ai, and c or 3 independent functions of them 
(e.g., cell volume V, shear angle $, and c/a). We define 
diagonal vector b = a 2 — ai, which is orthogonal to c 
and shuffle s (directed along the other diagonal ai +a 2 ), 


see Fig. [lj We emphasize that the complete description 
of this bcc-hcp transformation requires consideration of 
all 4 independent atomic and cell DoF. 

Previously, the bcc-hcp transformation path was usu¬ 
ally described by oversimplified models 0 0. In par¬ 
ticular, the shuffle-shear model [7 was a drag method 
that ignored critical changes in magnetization and vol¬ 
ume and hence improperly uncoupled atomic and cell 
DoF. A drag method improperly allows for discontinu¬ 
ity in some of the uncontrolled DoF; hence, the result¬ 
ing transformation path can “tunnel” under the TS (due 
to discontinuous DoF) and bypass it, missing and hence 
underestimating the true enthalpy barrier [T. Use of a 
drag method based on the rapid-nuclear-motion (RNM) 
approximation leads to a discontinuity in atomic shuf¬ 
fle, giving a very low bcc-hcp barrier in the form of a 
cusp, with equal enthalpies of the bcc and hep phases at 
the calculated pressure of 13.1 GPa [8]. We emphasize 
that the shuffle-shear model |7], which is a drag method, 
controls of only 2 DoF (shuffle and shear), leaving an un¬ 
controlled unit cell volume beyond its consideration |9]. 

The iron bcc-hcp transformation path was calculated 
directly using density functional theory (DFT) for the 
shuffle-shear model in a 2-atom cell (Fig. |T]) , with full 
volume relaxation for each pair of shuffle-shear values 
[9], finding a MEP with a cusp due to the magneto¬ 
volume collapse. In addition to a discontinuity in vol¬ 
ume, the reported MEP cusp was rounded due to the use 
of an approximate symmetry-adapted polynomial to an¬ 
alytically determine the entire potential energy surface. 
More recently, possible concerted transformations were 
investigated (including 3 shuffling mechanisms within the 
RNM at constant shear and fixed 71.5 bohr 3 /atom vol¬ 
ume) , finding a TS in the form of a cusp EH, but missing 
the observed magneto-volume collapse - disappearance of 
magnetization with a concomitant 11% drop in volume. 

Classical inter-atomic potentials typically ignore de¬ 
pendence on the magnetic state, giving an irrelevant path 
for magneto-structural transitions. Thus, the embedded 
atom method (EAM) provided an overestimated transi¬ 
tion pressure of 31-33 GPa for the uniform m and 14 
GPa for the uniaxial compression [TS .. 

While the pressure-induced bcc-hcp transformation in 
iron has been extensively studied, its complete theoret¬ 
ical description accounting for all relevant DoF is still 
needed, and we provide it below. Recall that the tra¬ 
ditional SSNEB method [4] can search the entire phase 
space, but requires a continuous and differentiable PES 
with the TS in the form of a saddle (e.g., as in Li [29]). 
which can be addressed by the available codes m with 
one m or (more stable) two [28] climbing images, while 
for a magneto-structural transformation the TS is a cusp 
on the PES with a discontinuity of atomic forces and 
stress components, where a climbing image will not stop. 
We extend the unrestricted SSNEB formalism (Appendix 
|A| ) and apply it to the bcc-hcp transformation in iron un¬ 
der pressure. We compare our results to experiment and 
contrast them with the previous theoretical discussions. 




















3 



P, GPa 

FIG. 4: (color online) Versus external pressure P, enthalpy 
[meV/Fe] of the TS (156 meV/Fe at Po =8.42 GPa) relative to 
bee (blue) and hep (red); unit cell volume of the TS (black), 
bcc (blue), and hep (red) structures; and shuffle at the TS. 


II. RESULTS 

Using our modified SSNEB method (see Appendix |A|) 
with properly coupled cell and atomic DoF, we perform 
calculations at several values of the applied hydrostatic 
pressure P, including 0, 8.4, 10, and 20 GPa. Because 
an equilibrium coexistence is rarely observed in experi¬ 
ment, here we provide the barrier versus pressure (Figs.[2j 
[3J and [4|, emphasizing generality of our quantitative 
description. The result at the equilibrium coexistence 
Po=8.4GPa is given in Fig. 3 in [6] (see also our Fig. [5]). 

A. Atomic and cell degrees of freedom 

The 2-atom cell has 6 atomic DoF (3 of which are in¬ 
dependent) and 6 cell DoF (3 lengths of the lattice trans¬ 
lation vectors and 3 angles between them). We find that 
out of the 3 independent atomic DoF, only 1 (shuffle) is 
interesting, while the other two can be chosen in such a 
way that they remain constant during the whole trans¬ 
formation. Among the 6 cell DoF, values of 2 out of 3 
angles remain constant (90°). The cell can be chosen with 
a mirror symmetry, so that 2 lengths of the translation 
vectors are equal (Fig. [I]). We are left with 1 changing 
atomic DoF (shuffle) and 3 cell DoF. The direct coordi¬ 
nates of atoms are (0, 0, 0) and (^[1 + §]; |[1 + §]; §), 
so that the second ones change from (1/2, 1/2, 1/2) in 
bcc to (2/3, 2/3, 1/2) in hep (Fig. [l]). 

We find that the atomic shuffle s is a continuous mono¬ 
tonic function of the SSNEB path variable, hence all the 
other path-dependent variables can be expressed in terms 
of the shuffle. All 3 independent cell DoF are plotted 
in terms of the shuffle in Fig. [3J where we introduced 



FIG. 5: (color online) For 2-atom cell, SSNEB enthalpy H 
[meV/atom], magnetization M [/is /atom], volume V [A 3 per 
cell], and internal pressure P [GPa], versus MEP. SSNEB re¬ 
sults for variable M (analytic continuation of H near cusp) 
are compared to those at fixed M, i.e., M — 0 (blue short- 
dash line) and 2.13/xs/Fe (red long-dash line). Two fixed-M 
solutions never cross on the PES, so there is a jump between 
them (vertical brown dashed line, as seen in V vs. path) to 
avoid higher-enthalpy solutions (blue and red dotted lines). 
Higher-enthalpy intermediate state (black squares) with in¬ 
termediate M arises from improper VASP convergence. 


b = 2a sin ($/2) as the length of b = a 2 — ai (Fig. [I]). 
Figures [2] and [3] show that although there is a discon¬ 
tinuity in magnetization at the TS, all the atomic and 
cell DoF are continuous functions versus the MEP (and 
the shuffle). This is in contrast to drag methods that 
have unphysical discontinuities in some of the uncon¬ 
trolled DoF. At Po = 8.4 GPa, the TS is characterized 
by a = 2.41 A, b = 2.59 A, c = 3.99 A, and s = 0.6, with 
V = 21.0 A 3 , 1? = 65.1°, and c/a = 1.6594. 

Although the SSNEB path depends on the Jacobian J, 
coupling atomic and cell degrees of freedom [4 , the TS is 
invariant for reasonable finite values of J (typically be¬ 
tween 1 and 5). At any J, the SSNEB path goes through 
the TS and the terminal states. With J & 3, we find 
an approximately linear dependence between the atomic 
shuffle s and cell translation b (Fig. |3|. Dependence of 
the shuffle at the TS on the hydrostatic external pressure 
P is nearly linear (Fig. [4] bottom). At the equilibrium 
coexistence pressure Po=8.42 GPa, where bcc and hep 
phases have equal enthalpies [6], the calculated enthalpy 
of the TS is 156 meV/atom above the terminal states. 

We emphasize that this value is the same from the in¬ 
tersection of lines connecting the SSNEB results at 0, 10, 
and 20 GPa (Fig. [4|, from our SSNEB calculation at P 0 
(Fig. |5|, and from direct DFT calculations of the elec¬ 
tronic structure (Fig. [6| and enthalpy of the TS atomic 
configuration at Po (Fig. [7]). 
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FIG. 6: (color online) Total spin DOS [states per atom per 
eV] for the bcc-hcp transformation at Po = 8.4 GPa. Major¬ 
ity (minority) spin states are plotted on positive (negative) 
axes. Endpoint are (a) non-magnetic hep and (e) magnetic 
(M = 2.12 /zb/F e) bcc. Transition states (TS) are (b) non¬ 
magnetic (M = 0) and (d) magnetic (M = 2 fis /Fe), having 
same atomic structure and enthalpy (the cusp in Fig. §. (c) 
Higher-enthalpy intermediate state (M=l.l /jlb/ Fe) with the 
same atomic structure as the TS. States (b, c, d) are en¬ 
thalpy extrema, marked (red squares) in Fig. [7] The right 
panels show interatomic bonds up to 2.5 A and 0.02 e~/A 3 
iso-surfaces of the spin density for the magnetic structures. 


B. Electronic structure and magnetization 

A magneto-structural transformation is characterized 
by a duality of the TS, with the same values of all the 
cell and atomic coordinates and the same enthalpy, but 
different magnetization (Fig. [2| and different electronic 
structure and density of states (DOS, Fig. [6|. Each in¬ 
termediate (higher enthalpy) state with a fractional mag¬ 
netization (filled shapes in Fig. [2]) can be regarded as an 
electronically excited state (Figf"[7]). The DOS of the in¬ 
termediate state is compared with DOS of two TS and 
two endpoints in Fig. [6j The bcc and hep end points 
and both TS have local minimum of DOS at the Fermi 
level. Bands continuously change between the end point 
and the corresponding TS, but the Fermi surface remains 
at the local minimum of the total DOS. Although the 
bands are not rigid, and the rigid-band model [32] is not 
precise here, the difference between two TS is roughly a 
band shift, with a local maximum in majority and minor¬ 
ity spins passing the Fermi surface. The higher-enthalpy 
intermediate state has that band at the Fermi surface. 

Magnetization is not a conserved DoF and is allowed 



FIG. 7: (color online) Enthalpy (meV/atom) vs. magnetiza¬ 
tion M (/ib/ atom) relative to the equal enthalpies of bcc and 
hep at Po = 8.4 GPa for the transition and intermediate states 
with the same atomic structure. The extrema are marked by 
red squares. The inset shows atomic structure with 0.2 e - /A 
iso-surfaces of the total electron density for the magnetic TS. 


to have discontinuity at the TS. DFT is expected to con¬ 
verge to a local enthalpy minimum for each atomic struc¬ 
ture, including the TS. Figure [7] shows two such minima 
(non-magnetic at M = 0 and magnetic at M^2/i#/Fe) 
with the same enthalpy but different values of magneti¬ 
zation; they are separated by the higher-enthalpy inter¬ 
mediate states, which can be regarded as electronically 
excited states with different values of M. Recall that 
both TS and intermediate states have the same atomic 
structure (identical atomic positions and lattice vectors), 
but different electronic structure. 

Notably, the directly calculated F vs. M (Fig. [7]) in¬ 
dicates that the states with intermediate magnetization 
(obtained by DFT with fixed M) have unstable electronic 
structure, but reflect an accommodation for the allowed 
discontinuity in M at the TS. Two electronically stable 
states at the H minima in Fig. [7] provide the enthalpy 
barrier, while all the others (including the maximal H ) 
do not, because they are electronically unstable. DFT 
is expected to find a local enthalpy minimum for any 
structure, including the TS; however, we find that the 
widely used DFT code VASP [33] [34] converges the elec¬ 
tronic structure to a local enthalpy extremum, which is 
maximum near the TS and minimum elsewhere. Careful 
application of the SSNEB approaching the TS (cusp re¬ 
gion) avoids those “intermediate” electronically excited 
states (filled shapes in Fig. [2| d ue to the additional stop 
condition, discussed in Sec. S and | A 3| 


C. Comparison to drag methods 

Most previous attempts to model the bcc-hcp transfor¬ 
mation in iron ultimately implemented a drag method, 
with discontinuity in one or more DoF. The RNM ap¬ 
proximation |8i results in a discontinuity in atomic shuf¬ 
fle. The shuffle-shear model [7] allows discontinuity of 
volume and lattice parameters [9]. By selecting too small 
(zero) or too large (infinite) Jacobian J within the SS¬ 
NEB [4], we can ignore either atomic (as in RNM) or 
cell DoF and reproduce the results of those drag meth- 
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ods with discontinuity in the ignored DoF, included in 
SSNEB with a negligible weight J-> 0 or 1/J< 1. 

With a discontinuity in any DoF, the physical system 
is dragged under an energy barrier; the true barrier is 
bypassed, and the calculated fictitious barrier is typically 
lower than the actual physical one. On the other hand, 
fixing or limiting selected DoF imposes restrictions on 
possible transition paths, constraining the search space, 
so that the minimum energy path can be missed: in this 
case the calculated barrier can be higher than the TS on 
the MEP. 

To avoid dealing with a cusp, it is tempting to con¬ 
sider two solutions with fixed magnetization (Fig.j5|. Al¬ 
though two projected energy curves versus path have ap¬ 
parent intersection in Fig. [5j this is not a TS, because two 
sets of DoF do not coincide anywhere. In other words, we 
find that fixed-magnetization solutions with M = 0 and 
2.13/ig/Fe do not intersect in a multi-dimensional space 
parametrized by all DoF. In particular, unit cell volume 
of the magnetic solution is larger than that of the non¬ 
magnetic one at every point in Fig. [5] The apparent 
intersection of the energy vs. path projected curves is a 
fictitious barrier (with an improper discontinuity in vol¬ 
ume and lattice constants), which is lower than the real 
one. At fixed magnetization, energy is a smooth function 
of DoF; and pressure and atomic forces can be converged 
to zero everywhere along such paths. 


D. Pressure discontinuity at the TS 

In contrast, pressure can be discontinuous at the TS. 
Two states with two different magnetization values are 
forced to have the same atomic and cell coordinates. 
Both are strained, and those strains have equal ampli¬ 
tudes and opposite directions in the dual TS. The TS 
relaxation to zero pressure would create discontinuity in 
one or more DoF (this was previously done in drag meth¬ 
ods). In particular, continuity of volume creates discon¬ 
tinuity of pressure, and vice versus. 

Interestingly, pressure discontinuity of 25 GPa (from 
— 13 to 12 GPa in Fig. [5| at the TS has the same or¬ 
der of magnitude as the spread of experimental pressures 
for coexisting bcc and hep phases (from 0 to 25 GPa). 
This is not a coincidence. Additional broadening of the 
measured pressures is caused by the martensitic stress. 

Indeed, the bcc-hcp (as) transition in iron is marten¬ 
sitic, and the hysteresis loop can be characterized by four 
pressure values: Pgta^t^ an d ^end^ f° r direct bcc->hcp; 

Pstart ^ an d Pend ^ f° r the reverse hep—)Ucc transform 
[6]. In one experiment (7], the hep phase appears at 

Pstart^ = 10-8 ± 0.5 GPa and bcc disappears above 
= 21 GPa upon loading; upon unloading, bcc 
reappears at Pgt^t* = 15.8 ±0.5 GPa and hep disappears 
below P^2°^ = 3 GPa . The observed bcc-to-hcp onset 
Pstart'* das been reported from 8.6 to 15 GPa, with the 


highest experimental P^2°^ = 8.5 ± 0.6 GPa 0. Our 
calculated hydrostatic equilibrium coexistence pressure 
of 8.4 GPa agrees well with a range of experimental mea¬ 
surements [6]. The 25 GPa pressure discontinuity at the 
TS agrees with the experimental « 25 GPa spread of the 
observed bcc-hcp coexistence (Table 1 in [6]). 


III. SUMMARY 


In summary, we extended the SSNEB formalism to ad¬ 
dress magneto-volume collapse in magnetic materials un¬ 
der stress or temperature, caused by loss of magnetiza¬ 
tion along the transformation pathway, where magnetiza¬ 
tion is a non-conserved order parameter, unlike the con¬ 
tinuous atomic and cell degrees of freedom. We applied 
it to the pressure-induced bcc-hcp magneto-structural 
transformation in iron, which exhibits a well-known 11% 
volume decrease concomitant with the loss of magneti¬ 
zation. We contrasted DFT-based equilibrium coexis¬ 
tence pressure (typically found by Maxwell construction) 
and “apparent coexistence” pressures typically measured, 
which are highly affected by internal stresses in the sam¬ 
ples. For iron under pressure, we found the transition 
state in the form of a cusp on the potential energy sur¬ 
face, although all the atomic and cell degrees of freedom 
are continuous, as physically required. We explained the 
difference between the generalized SSNEB and all previ¬ 
ously used drag and restricted parametrization methods. 

Our calculated values of the bcc-hcp equilibrium co¬ 
existence pressure, energy barrier, and pressure discon¬ 
tinuity at the transition state all agree with the exper¬ 
imental data. The new formalism is suitable for study¬ 
ing numerous magnetic systems, especially those involv¬ 
ing magneto-structural transformations, such as trans¬ 
ducers, magnetic switches, and competing chemical and 
magnetic structures. Our extended formalism works for 
barriers in the forms of both saddle points and cusps. 
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Appendix A: Computational Methods 


We combine density functional theory (DFTjA 1|) wit h 
the extended nudged elastic band methods ( |A 2 


and 


A3). 
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1. DFT 


The Vienna ab initio simulation package VASP [33] [34] 
is used to calculate electronic energy, pressure, and 
atomic forces for instantaneous atomic configurations. 
We use the projector augmented wave (PAW) [H] tech¬ 
nique, the generalized gradient approximation (GGA 
m), a plane-wave energy cutoff of 700 eV, and a con¬ 
verged fc-point mesh [377 including T-point for the Bril- 
louin zone integration with at least 75 points per A -1 
(e.g., 17 ^-points for b = 4.4A). The modified Broy- 
den method [38] is used for self-consistency. Atomic 
and cell relaxations (including those within the NEB 
method) are performed by selective dynamics using a 
conjugate-gradient algorithm with a Gaussian smearing 
with a = 0.05 eV. Structural energy differences are ob¬ 
tained using the tetrahedron method with Blochl correc¬ 
tions. 

For any atomic structure, DFT should find a local en¬ 
ergy minimum in terms of the electronic density. VASP 
does it for most structures, except those near the TS. At 
the TS it converges to a local enthalpy maximum (Fig. [ 7 ]) 
- an intermediate state with a higher enthalpy and a 
fractional magnetization (filled symbols in Figs. [2] and 
[5|. Convergence to the nearest extremum (which is ex¬ 
pected to be a minimum) provides an unexpected result 
(maximum) near the TS. At present, we have no expla¬ 
nation for this feature. 


2. NEB with recursively added images near TS 

We recursively add images near the TS within the SS- 
NEB method [4]. First, a traditional SSNEB [4] with 
equidistant images is used. The input configuration is 
prepared by adjusting initial unit cell volume of each im¬ 
age in such a way, that all images between the magnetic 
bee end point and the expected TS are ferromagnetic 
(FM), while all images on the other side (TS to hep) 
have zero magnetization. After convergence of this SS¬ 
NEB calculation, we take two closest to the TS neighbor 
images with different magnetic states, and use them as 
new terminal points for the next SSNEB calculation with 
a sufficiently large number of images. Adding images be¬ 
tween those closest to the TS is repeated until the cusp 
is well-sampled by a dense grid. We stop these iterations 
when any image added between the two TS neighbors 
(one is ferromagnetic, the other is non-magnetic) gets an 
intermediate magnetization 0.1 < M < 1.8/iB/atom. 


3. NEB with two images approaching the TS 

The traditional SSNEB method [4] expects a TS to be 
a saddle point of the smooth potential energy surface. In 
the single climbing image algorithm [31], the highest en¬ 


thalpy image is driven up to the saddle point. This image 
does not feel the spring forces along the band. It tries to 
maximize its energy along the NEB, and minimize it in 
all other directions. This algorithm requires continuity 
of forces and pressures along the path, and expects zero 
forces at the saddle point; it is not suitable for PES with 
a barrier in the form of a cusp, where discontinuity of 
pressure and atomic forces is at issue. 

In contrast, magnetic and magneto-structural transfor¬ 
mations can be addressed by a modified SSNEB, where 
the TS is approached by two images with different mag¬ 
netization (Fig. 2|, converging towards the same values 
of all atomic and cell DoF and the same enthalpy. Ferro¬ 
magnetic and non-magnetic states form two smooth en¬ 
ergy surfaces; the minimal energy surface forms a cusp at 
their intersection; components of the atomic forces and 
stress tensor have discontinuities at that cusp (Fig. [5|. 

In our algorithm, except for the two terminal images 
and two images near the TS, all the other images are 
nudged. At each step, a trial image is placed between 
the two images with true tangent force projections hav¬ 
ing opposite directions (those two images bracket the TS 
from two sides). Initial coordinates of the trial image are 
the weighted average of coordinates of those two images: 
1-Trial = m/|+(1— m)l 2 , where the weight 0 < m < 1 can 
have any value between zero and one (we chose m = 1/2). 
The trial image is relaxed to the point with zero (or small, 
below the specified cutoff) perpendicular forces. The true 
tangent force projection is calculated for the trial image 
and its direction is compared with those for its two neigh¬ 
bors. The image with tangent force in the same direction 
is replaced by the trial image, so that the tangent forces 
of two climbing images have opposite directions again. 
Calculation stops if magnetization of any image has an 
unacceptable intermediate value, or the “distance” be¬ 
tween two images is below a certain cutoff. The normal¬ 
ized “distance” between images is the square root of the 
sum of weighted squares of the differences between com¬ 
ponents of atomic coordinates, lattice translation vectors, 
and two enthalpies. At each step, only one trial image is 
computed and only one of the two TS neighbors is moved. 

Our algorithm works for barriers in the forms of both 
saddle points and cusps. It can safely replace the 
climbing-image algorithm C 1-NEB |3T in all cases, where 
C 1-NEB is applicable, and can be generalized for complex 
energy landscapes [28]. However, the SSNEB conver¬ 
gence malfunctions if any image gets an excited electronic 
state with intermediate magnetization (Fig. [6| . Fortu¬ 
nately, such cases can be easily detected, and we monitor 
them by an additional check: the calculation stops if any 
image has an unacceptable intermediate magnetization, 
0.1 < M < 1.8/i#/citom. As a result, we get two im¬ 
ages with very similar values of enthalpy and all atomic 
and cell DoF, and an intermediate state with a higher 
enthalpy between them. A similar result was obtained 


from a dense grid of images near the TS in section A 2 




7 


[1] S. V. Vonsovsky, Magnetism (Nauka, Moscow, 1971). 

[2] A. Akhiezer, General physics. The electrical and mag¬ 
netic phenomena (Naukova Dumka, Kiev, 1981). 

[3] H. Jonsson, G. Mills and K. W. Jacobsen, Nudged Elas¬ 
tic Band Method for Finding Minimum Energy Paths of 
Transitions (World Scientific, 1998). 

[4] D. Sheppard, P. H. Xiao, W. Chemelewski, D. D. John¬ 
son, and G. Henkelman, J. Chem. Phys. 136, 074103 
( 2012 ). 

[5] V. V. Cherdyntsev, Abundance of Chemical Elements 
(University of Chicago Press, Chicago, 1961), translated 
from Russian (GosTechlzdat, Moscow, 1956). 

[6] N. A. Zarkevich and D. D. Johnson, Phys. Rev. B 91, 
174104 (2015). 

[7] W. A. Bassett and E. Huang, Science 238, 780 (1987). 

[8] D. F. Johnson and E. A. Carter, J. Chem. Phys. 128, 
104703 (2008). 

[9] J. B. Liu and D. D. Johnson, Phys. Rev. B 79, 134113 
(2009). 

[10] R. D. Taylor, M. P. Pasternak, and R. Jeanloz, Journal 
of Applied Physics 69, 6126 (1991). 

[11] J. Jung, M. Fricke, G. Hampel, and J. Hesse, Hyperhne 
Interactions 72, 375 (1992). 

[12] M. Ekman, B. Sadigh, K. Einarsdotter, and P. Blaha, 
Phys. Rev. B 58, 5296 (1998). 

[13] F. M. Wang and R. Ingalls, Phys. Rev. B 57, 5647 (1998). 

[14] M. Sob, M. Friak, L. G. Wang, and V. Vitek, Multiscale 
Modelling of Materials 538, 523 (1999). 

[15] M. Friak and M. Sob, Phys. Rev. B 77, 174117 (2008). 

[16] J. Z. Jiang, J. S. Olsen, and L. Gerward, Materials Trans¬ 
actions 42, 1571 (2001). 

[17] J. L. Shao, A. M. He, S. Q. Duan, P. Wang, and C. S. 
Qin, Acta Physica Sinica 59, 4888 (2010). 

[18] B. T. Wang, J. L. Shao, G. C. Zhang, W. D. Li, and 
P. Zhang, J. Phys.: Condens. Matter 21, 495702 (2009). 

[19] B. T. Wang, J. L. Shao, G. C. Zhang, W. D. Li, and 
P. Zhang, J. Phys.: Condens. Matter 22, 435404 (2010). 

[20] M. Ortiz, Abstracts of Papers of the American Chemical 


Society 233, 118 (2007). 

[21] B. Yaakobi, T. R. Boehly, D. D. Meyerhofer, T. J. B. 
Collins, B. A. Remington, P. G. Allen, S. M. Pollaine, 
H. E. Lorenzana, and J. H. Eggert, Phys. Rev. Lett. 95, 
075501 (2005). 

[22] F. Baudelet, S. Pascarelli, O. Mathon, J. P. Itie, A. Po- 
lian, M. d’ Astuto, and J. C. Chervin, J. Phys.: Condens. 
Matter 17, S957 (2005). 

[23] K. J. Caspersen, A. Lew, M. Ortiz, and E. A. Carter, 
Physical Review Letters 93, 115501 (2004). 

[24] N. Suzuki, T. Souraku, I. Hamada, and T. Takezawa, 
High Pressure Research 22, 451 (2002). 

[25] L. Stixrude and R. E. Cohen, Science 267, 1972 (1995). 

[26] L. Stixrude, E. Wasserman, and R. E. Cohen, Properties 
of Earth and Planetary Materials at High Pressure and 
Temperature 101, 159 (1998). 

[27] B. Dupe, B. Amadon, Y. P. Pellegrini, and C. Denoual, 
Phys. Rev. B 87, 024103 (2013). 

[28] N. A. Zarkevich and D. D. Johnson, J. Chem. Phys. 142, 
024106 (2015). 

[29] K. J. Caspersen and E. A. Carter, PNAS 102, 6738 
(2005). 

[30] N. A. Zarkevich and D. D. Johnson, C2NEB source code 
(2014), URL http://lib.dr.iastate.edu/ameslab_ 
software/1/ 

[31] G. Henkelman, B. P. Uberuaga, and H. Jonsson, J. Chem. 
Phys. 113, 9901 (2000). 

[32] N. F. Mott and H. Jones, The Theory of the Properties 
of Metals and Alloys (Dover Publications, New York,, 
1958). 

[33] G. Kresse and J. Hafner, Phys. Rev. B 47, RC558 (1993). 

[34] G. Kresse, and J. Hafner, Phys. Rev. B 49, 14251 (1994). 

[35] P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

[36] J. P. Perdew, Phys. Lett. A 165, 79 (1992). 

[37] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 
(1976). 

[38] D. D. Johnson, Phys. Rev. B 38, 12807 (1988). 



